home *** CD-ROM | disk | FTP | other *** search
/ Amiga Format CD 33 / Amiga Format AFCD33 (Issue 117, Dec 1998).iso / -seriously_amiga- / sound / mpeginoutppc / tonal.c < prev   
C/C++ Source or Header  |  1998-09-07  |  37KB  |  1,035 lines
  1. /**********************************************************************
  2. Copyright (c) 1991 MPEG/audio software simulation group, All Rights Reserved
  3. tonal.c
  4. **********************************************************************/
  5. /**********************************************************************
  6.  * MPEG/audio coding/decoding software, work in progress              *
  7.  *   NOT for public distribution until verified and approved by the   *
  8.  *   MPEG/audio committee.  For further information, please contact   *
  9.  *   Davis Pan, 508-493-2241, e-mail: pan@3d.enet.dec.com             *
  10.  *                                                                    *
  11.  * VERSION 3.9t                                                       *
  12.  *   changes made since last update:                                  *
  13.  *   date   programmers         comment                               *
  14.  * 2/25/91  Douglas Wong        start of version 1.1 records          *
  15.  * 3/06/91  Douglas Wong        rename: setup.h to endef.h            *
  16.  *                              updated I_psycho_one and II_psycho_one*
  17.  * 3/11/91  W. J. Carter        Added Douglas Wong's updates dated    *
  18.  *                              3/9/91 for I_Psycho_One() and for     *
  19.  *                              II_Psycho_One().                      *
  20.  * 5/10/91  W. Joseph Carter    Ported to Macintosh and Unix.         *
  21.  *                              Located and fixed numerous software   *
  22.  *                              bugs and table data errors.           *
  23.  * 6/11/91  Davis Pan           corrected several bugs                *
  24.  *                              based on comments from H. Fuchs       *
  25.  * 01jul91  dpwe (Aware Inc.)   Made pow() args float                 *
  26.  *                              Removed logical bug in I_tonal_label: *
  27.  *                              Sometimes *tone returned == STOP      *
  28.  * 7/10/91  Earle Jennings      no change necessary in port to MsDos  *
  29.  * 11sep91  dpwe@aware.com      Subtracted 90.3dB from II_f_f_t peaks *
  30.  * 10/1/91  Peter W. Farrett    Updated II_Psycho_One(),I_Psycho_One()*
  31.  *                              to include comments.                  *
  32.  *11/29/91  Masahiro Iwadare    Bug fix regarding POWERNORM           *
  33.  *                              fixed several other miscellaneous bugs*
  34.  * 2/11/92  W. Joseph Carter    Ported new code to Macintosh.  Most   *
  35.  *                              important fixes involved changing     *
  36.  *                              16-bit ints to long or unsigned in    *
  37.  *                              bit alloc routines for quant of 65535 *
  38.  *                              and passing proper function args.     *
  39.  *                              Removed "Other Joint Stereo" option   *
  40.  *                              and made bitrate be total channel     *
  41.  *                              bitrate, irrespective of the mode.    *
  42.  *                              Fixed many small bugs & reorganized.  *
  43.  * 2/12/92  Masahiro Iwadare    Fixed some potential bugs in          *
  44.  *          Davis Pan           subsampling()                         *
  45.  * 2/25/92  Masahiro Iwadare    Fixed some more potential bugs        *
  46.  * 6/24/92  Tan Ah Peng         Modified window for FFT               * 
  47.  *                              (denominator N-1 to N)                *
  48.  *                              Updated all critical band rate &      *
  49.  *                              absolute threshold tables and critical*
  50.  *                              boundaries for use with Layer I & II  *  
  51.  *                              Corrected boundary limits for tonal   *
  52.  *                              component computation                 *
  53.  *                              Placement of non-tonal component at   *
  54.  *                              geometric mean of critical band       *
  55.  *                              (previous placement method commented  *
  56.  *                               out - can be used if desired)        *
  57.  * 3/01/93  Mike Li             Infinite looping fix in noise_label() *
  58.  * 3/19/93  Jens Spille         fixed integer overflow problem in     *
  59.  *                              psychoacoutic model 1                 *
  60.  * 3/19/93  Giorgio Dimino      modifications to better account for   *
  61.  *                              tonal and non-tonal components        *
  62.  * 5/28/93 Sriram Jayasimha     "London" mod. to psychoacoustic model1*
  63.  * 8/05/93 Masahiro Iwadare     noise_label modification "option"     *
  64.  **********************************************************************/
  65.  
  66. #include "common.h"
  67. #include "encoder.h"
  68. #define LONDON                  /* enable "LONDON" modification */
  69. #define MAKE_SENSE              /* enable "MAKE_SENSE" modification */
  70. #define MI_OPTION               /* enable "MI_OPTION" modification */
  71. /**********************************************************************/
  72. /*
  73. /*        This module implements the psychoacoustic model I for the
  74. /* MPEG encoder layer II. It uses simplified tonal and noise masking
  75. /* threshold analysis to generate SMR for the encoder bit allocation
  76. /* routine.
  77. /*
  78. /**********************************************************************/
  79.  
  80. int crit_band;
  81. int FAR *cbound;
  82. int sub_size;
  83.  
  84. void read_cbound(lay,freq)  /* this function reads in critical */
  85. int lay, freq;              /* band boundaries                 */
  86. {
  87.  int i,j,k;
  88.  FILE *fp;
  89.  char r[16], t[80];
  90.  
  91.  strcpy(r, "2cb1");
  92.  r[0] = (char) lay + '0';
  93.  r[3] = (char) freq + '0';
  94.  if( !(fp = OpenTableFile(r)) ){       /* check boundary values */
  95.     printf("Please check %s boundary table\n",r);
  96.     exit(1);
  97.  }
  98.  fgets(t,80,fp);               /* read input for critical bands */
  99.  sscanf(t,"%d\n",&crit_band);
  100.  cbound = (int FAR *) mem_alloc(sizeof(int) * crit_band, "cbound");
  101.  for(i=0;i<crit_band;i++){   /* continue to read input for */
  102.     fgets(t,80,fp);            /* critical band boundaries   */
  103.     sscanf(t,"%d %d\n",&j, &k);
  104.     if(i==j) cbound[j] = k;
  105.     else {                     /* error */
  106.        printf("Please check index %d in cbound table %s\n",i,r);
  107.        exit(1);
  108.     }
  109.  }
  110.  fclose(fp);
  111. }        
  112.  
  113. void read_freq_band(ltg,lay,freq)  /* this function reads in   */
  114. int lay, freq;                     /* frequency bands and bark */
  115. g_ptr FAR *ltg;                /* values                   */
  116. {
  117.  int i,j, k;
  118.  double b,c;
  119.  FILE *fp;
  120.  char r[16], t[80];
  121.  
  122.  strcpy(r, "2th1");
  123.  r[0] = (char) lay + '0';
  124.  r[3] = (char) freq + '0';
  125.  if( !(fp = OpenTableFile(r)) ){   /* check freq. values  */
  126.     printf("Please check frequency and cband table %s\n",r);
  127.     exit(1);
  128.  }
  129.  fgets(t,80,fp);              /* read input for freq. subbands */
  130.  sscanf(t,"%d\n",&sub_size);
  131.  *ltg = (g_ptr FAR ) mem_alloc(sizeof(g_thres) * sub_size, "ltg");
  132.  (*ltg)[0].line = 0;          /* initialize global masking threshold */
  133.  (*ltg)[0].bark = 0;
  134.  (*ltg)[0].hear = 0;
  135.  for(i=1;i<sub_size;i++){    /* continue to read freq. subband */
  136.     fgets(t,80,fp);          /* and assign                     */
  137.     sscanf(t,"%d %d %lf %lf\n",&j, &k, &b, &c);
  138.     if(i == j){
  139.        (*ltg)[j].line = k;
  140.        (*ltg)[j].bark = b;
  141.        (*ltg)[j].hear = c;
  142.     }
  143.     else {                   /* error */
  144.        printf("Please check index %d in freq-cb table %s\n",i,r);
  145.        exit(1);
  146.     }
  147.  }
  148.  fclose(fp);
  149. }
  150.  
  151. void make_map(power, ltg)       /* this function calculates the */
  152. mask FAR power[HAN_SIZE];   /* global masking threshold     */
  153. g_thres FAR *ltg;
  154. {
  155.  int i,j;
  156.  
  157.  for(i=1;i<sub_size;i++) for(j=ltg[i-1].line;j<=ltg[i].line;j++)
  158.     power[j].map = i;
  159. }
  160.  
  161. double add_db(a,b)
  162. double a,b;
  163. {
  164.  a = pow(10.0,a/10.0);
  165.  b = pow(10.0,b/10.0);
  166.  return 10 * log10(a+b);
  167. }
  168.  
  169. /****************************************************************/
  170. /*
  171. /*        Fast Fourier transform of the input samples.
  172. /*
  173. /****************************************************************/
  174.  
  175. void II_f_f_t(sample, power)      /* this function calculates an */
  176. double FAR sample[FFT_SIZE];  /* FFT analysis for the freq.  */
  177. mask FAR power[HAN_SIZE];     /* domain                      */
  178. {
  179.  int i,j,k,L,l=0;
  180.  int ip, le, le1;
  181.  double t_r, t_i, u_r, u_i;
  182.  static int M, MM1, init = 0, N;
  183.  double *x_r, *x_i, *energy;
  184.  static int *rev;
  185.  static double *w_r, *w_i;
  186.  
  187.  x_r = (double *) mem_alloc(sizeof(DFFT), "x_r");
  188.  x_i = (double *) mem_alloc(sizeof(DFFT), "x_i");
  189.  energy = (double *) mem_alloc(sizeof(DFFT), "energy");
  190.  for(i=0;i<FFT_SIZE;i++) x_r[i] = x_i[i] = energy[i] = 0;
  191.  if(!init){
  192.     rev = (int *) mem_alloc(sizeof(IFFT), "rev");
  193.     w_r = (double *) mem_alloc(sizeof(D10), "w_r");
  194.     w_i = (double *) mem_alloc(sizeof(D10), "w_i");
  195.     M = 10;
  196.     MM1 = 9;
  197.     N = FFT_SIZE;
  198.     for(L=0;L<M;L++){
  199.        le = 1 << (M-L);
  200.        le1 = le >> 1;
  201.        w_r[L] = cos(PI/le1);
  202.        w_i[L] = -sin(PI/le1);
  203.     }
  204.     for(i=0;i<FFT_SIZE;rev[i] = l,i++) for(j=0,l=0;j<10;j++){
  205.        k=(i>>j) & 1;
  206.        l |= (k<<9-j);                
  207.     }
  208.     init = 1;
  209.  }
  210.  memcpy( (char *) x_r, (char *) sample, sizeof(double) * FFT_SIZE);
  211.  for(L=0;L<MM1;L++){
  212.     le = 1 << (M-L);
  213.     le1 = le >> 1;
  214.     u_r = 1;
  215.     u_i = 0;
  216.     for(j=0;j<le1;j++){
  217.        for(i=j;i<N;i+=le){
  218.           ip = i + le1;
  219.           t_r = x_r[i] + x_r[ip];
  220.           t_i = x_i[i] + x_i[ip];
  221.           x_r[ip] = x_r[i] - x_r[ip];
  222.           x_i[ip] = x_i[i] - x_i[ip];
  223.           x_r[i] = t_r;
  224.           x_i[i] = t_i;
  225.           t_r = x_r[ip];
  226.           x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
  227.           x_i[ip] = x_i[ip] * u_r + t_r * u_i;
  228.        }
  229.        t_r = u_r;
  230.        u_r = u_r * w_r[L] - u_i * w_i[L];
  231.        u_i = u_i * w_r[L] + t_r * w_i[L];
  232.     }
  233.  }
  234.  for(i=0;i<N;i+=2){
  235.     ip = i + 1;
  236.     t_r = x_r[i] + x_r[ip];
  237.     t_i = x_i[i] + x_i[ip];
  238.     x_r[ip] = x_r[i] - x_r[ip];
  239.     x_i[ip] = x_i[i] - x_i[ip];
  240.     x_r[i] = t_r;
  241.     x_i[i] = t_i;
  242.     energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
  243.  }
  244.  for(i=0;i<FFT_SIZE;i++) if(i<rev[i]){
  245.     t_r = energy[i];
  246.     energy[i] = energy[rev[i]];
  247.     energy[rev[i]] = t_r;
  248.  }
  249.  for(i=0;i<HAN_SIZE;i++){    /* calculate power density spectrum */
  250.     if (energy[i] < 1E-20) energy[i] = 1E-20;
  251.     power[i].x = 10 * log10(energy[i]) + POWERNORM;
  252.     power[i].next = STOP;
  253.     power[i].type = FALSE;
  254.  }
  255.  mem_free((void **) &x_r);
  256.  mem_free((void **) &x_i);
  257.  mem_free((void **) &energy);
  258. }
  259.  
  260. /****************************************************************/
  261. /*
  262. /*         Window the incoming audio signal.
  263. /*
  264. /****************************************************************/
  265.  
  266. void II_hann_win(sample)          /* this function calculates a  */
  267. double FAR sample[FFT_SIZE];  /* Hann window for PCM (input) */
  268. {                                 /* samples for a 1024-pt. FFT  */
  269.  register int i;
  270.  register double sqrt_8_over_3;
  271.  static int init = 0;
  272.  static double FAR *window;
  273.  
  274.  if(!init){  /* calculate window function for the Fourier transform */
  275.     window = (double FAR *) mem_alloc(sizeof(DFFT), "window");
  276.     sqrt_8_over_3 = pow(8.0/3.0, 0.5);
  277.     for(i=0;i<FFT_SIZE;i++){
  278.        /* Hann window formula */
  279.        window[i]=sqrt_8_over_3*0.5*(1-cos(2.0*PI*i/(FFT_SIZE)))/FFT_SIZE;
  280.     }
  281.     init = 1;
  282.  }
  283.  for(i=0;i<FFT_SIZE;i++) sample[i] *= window[i];
  284. }
  285.  
  286. /*******************************************************************/
  287. /*
  288. /*        This function finds the maximum spectral component in each
  289. /* subband and return them to the encoder for time-domain threshold
  290. /* determination.
  291. /*
  292. /*******************************************************************/
  293. #ifndef LONDON
  294. void II_pick_max(power, spike)
  295. double FAR spike[SBLIMIT];
  296. mask FAR power[HAN_SIZE];
  297. {
  298.  double max;
  299.  int i,j;
  300.  
  301.  for(i=0;i<HAN_SIZE;spike[i>>4] = max, i+=16)      /* calculate the      */
  302.  for(j=0, max = DBMIN;j<16;j++)                    /* maximum spectral   */
  303.     max = (max>power[i+j].x) ? max : power[i+j].x; /* component in each  */
  304. }                                                  /* subband from bound */
  305.                                                    /* 4-16               */
  306. #else
  307. void II_pick_max(power, spike)
  308. double FAR spike[SBLIMIT];
  309. mask FAR power[HAN_SIZE];
  310. {
  311.  double sum;
  312.  int i,j;
  313.  
  314.  for(i=0;i<HAN_SIZE;spike[i>>4] = 10.0*log10(sum), i+=16)
  315.                                                    /* calculate the      */
  316.  for(j=0, sum = pow(10.0,0.1*DBMIN);j<16;j++)      /* sum of spectral   */
  317.    sum += pow(10.0,0.1*power[i+j].x);              /* component in each  */
  318. }                                                  /* subband from bound */
  319.                                                    /* 4-16               */
  320. #endif
  321.  
  322. /****************************************************************/
  323. /*
  324. /*        This function labels the tonal component in the power
  325. /* spectrum.
  326. /*
  327. /****************************************************************/
  328.  
  329. void II_tonal_label(power, tone)  /* this function extracts (tonal) */
  330. mask FAR power[HAN_SIZE];     /* sinusoidals from the spectrum  */
  331. int *tone;
  332. {
  333.  int i,j, last = LAST, first, run, last_but_one = LAST; /* dpwe */
  334.  double max;
  335.  
  336.  *tone = LAST;
  337.  for(i=2;i<HAN_SIZE-12;i++){
  338.     if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
  339.        power[i].type = TONE;
  340.        power[i].next = LAST;
  341.        if(last != LAST) power[last].next = i;
  342.        else first = *tone = i;
  343.        last = i;
  344.     }
  345.  }
  346.  last = LAST;
  347.  first = *tone;
  348.  *tone = LAST;
  349.  while(first != LAST){               /* the conditions for the tonal          */
  350.     if(first<3 || first>500) run = 0;/* otherwise k+/-j will be out of bounds */
  351.     else if(first<63) run = 2;       /* components in layer II, which         */
  352.     else if(first<127) run = 3;      /* are the boundaries for calc.          */
  353.     else if(first<255) run = 6;      /* the tonal components                  */
  354.     else run = 12;
  355.     max = power[first].x - 7;        /* after calculation of tonal   */
  356.     for(j=2;j<=run;j++)              /* components, set to local max */
  357.        if(max < power[first-j].x || max < power[first+j].x){
  358.           power[first].type = FALSE;
  359.           break;
  360.        }
  361.     if(power[first].type == TONE){   /* extract tonal components */
  362.        int help=first;
  363.        if(*tone==LAST) *tone = first;
  364.        while((power[help].next!=LAST)&&(power[help].next-first)<=run)
  365.           help=power[help].next;
  366.        help=power[help].next;
  367.        power[first].next=help;
  368.        if((first-last)<=run){
  369.           if(last_but_one != LAST) power[last_but_one].next=first;
  370.        }
  371.        if(first>1 && first<500){     /* calculate the sum of the */
  372.           double tmp;                /* powers of the components */
  373.           tmp = add_db(power[first-1].x, power[first+1].x);
  374.           power[first].x = add_db(power[first].x, tmp);
  375.        }
  376.        for(j=1;j<=run;j++){
  377.           power[first-j].x = power[first+j].x = DBMIN;
  378.           power[first-j].next = power[first+j].next = STOP;
  379.           power[first-j].type = power[first+j].type = FALSE;
  380.        }
  381.        last_but_one=last;
  382.        last = first;
  383.        first = power[first].next;
  384.     }
  385.     else {
  386.        int ll;
  387.        if(last == LAST); /* *tone = power[first].next; dpwe */
  388.        else power[last].next = power[first].next;
  389.        ll = first;
  390.        first = power[first].next;
  391.        power[ll].next = STOP;
  392.     }
  393.  }
  394. }
  395.  
  396. /****************************************************************/
  397. /*
  398. /*        This function groups all the remaining non-tonal
  399. /* spectral lines into critical band where they are replaced by
  400. /* one single line.
  401. /*
  402. /****************************************************************/
  403.         
  404. void noise_label(power, noise, ltg)
  405. g_thres FAR *ltg;
  406. mask FAR *power;
  407. int *noise;
  408. {
  409.  int i,j, centre, last = LAST;
  410.  double index, weight, sum;
  411.                               /* calculate the remaining spectral */
  412.  for(i=0;i<crit_band-1;i++){  /* lines for non-tonal components   */
  413.      for(j=cbound[i],weight = 0.0,sum = DBMIN;j<cbound[i+1];j++){
  414.         if(power[j].type != TONE){
  415.            if(power[j].x != DBMIN){
  416.               sum = add_db(power[j].x,sum);
  417. /* the line below and others under the "MAKE_SENSE" condition are an alternate
  418.    interpretation of "geometric mean". This approach may make more sense but
  419.    it has not been tested with hardware. */
  420. #ifdef MAKE_SENSE
  421.               weight += pow(10.0, power[j].x/10.0) * (ltg[power[j].map].bark-i);
  422. #endif
  423.               power[j].x = DBMIN;
  424.            }
  425.         }   /*  check to see if the spectral line is low dB, and if  */
  426.      }      /* so replace the center of the critical band, which is */
  427.             /* the center freq. of the noise component              */
  428.  
  429. #ifdef MAKE_SENSE
  430.      if(sum <= DBMIN)  centre = (cbound[i+1]+cbound[i]) /2;
  431.      else {
  432.         index = weight/pow(10.0,sum/10.0);
  433.         centre = cbound[i] + (int) (index * (double) (cbound[i+1]-cbound[i]) );
  434.      } 
  435. #else
  436.      index = (double)( ((double)cbound[i]) * ((double)(cbound[i+1]-1)) );
  437.      centre = (int)(pow(index,0.5)+0.5);
  438. #endif
  439.  
  440.     /* locate next non-tonal component until finished; */
  441.     /* add to list of non-tonal components             */
  442. #ifdef MI_OPTION
  443.      /* Masahiro Iwadare's fix for infinite looping problem? */
  444.      if(power[centre].type == TONE) 
  445.        if (power[centre+1].type == TONE) centre++; else centre--;
  446. #else
  447.      /* Mike Li's fix for infinite looping problem */
  448.      if(power[centre].type == FALSE) centre++;
  449.  
  450.      if(power[centre].type == NOISE){
  451.        if(power[centre].x >= ltg[power[i].map].hear){
  452.          if(sum >= ltg[power[i].map].hear) sum = add_db(power[j].x,sum);
  453.          else
  454.          sum = power[centre].x;
  455.        }
  456.      }
  457. #endif
  458.      if(last == LAST) *noise = centre;
  459.      else {
  460.         power[centre].next = LAST;
  461.         power[last].next = centre;
  462.      }
  463.      power[centre].x = sum;
  464.      power[centre].type = NOISE;        
  465.      last = centre;
  466.  }        
  467. }
  468.  
  469. /****************************************************************/
  470. /*
  471. /*        This function reduces the number of noise and tonal
  472. /* component for further threshold analysis.
  473. /*
  474. /****************************************************************/
  475.  
  476. void subsampling(power, ltg, tone, noise)
  477. mask FAR power[HAN_SIZE];
  478. g_thres FAR *ltg;
  479. int *tone, *noise;
  480. {
  481.  int i, old;
  482.  
  483.  i = *tone; old = STOP;    /* calculate tonal components for */
  484.  while(i!=LAST){           /* reduction of spectral lines    */
  485.     if(power[i].x < ltg[power[i].map].hear){
  486.        power[i].type = FALSE;
  487.        power[i].x = DBMIN;
  488.        if(old == STOP) *tone = power[i].next;
  489.        else power[old].next = power[i].next;
  490.     }
  491.     else old = i;
  492.     i = power[i].next;
  493.  }
  494.  i = *noise; old = STOP;    /* calculate non-tonal components for */
  495.  while(i!=LAST){            /* reduction of spectral lines        */
  496.     if(power[i].x < ltg[power[i].map].hear){
  497.        power[i].type = FALSE;
  498.        power[i].x = DBMIN;
  499.        if(old == STOP) *noise = power[i].next;
  500.        else power[old].next = power[i].next;
  501.     }
  502.     else old = i;
  503.     i = power[i].next;
  504.  }
  505.  i = *tone; old = STOP;
  506.  while(i != LAST){                              /* if more than one */
  507.     if(power[i].next == LAST)break;             /* tonal component  */
  508.     if(ltg[power[power[i].next].map].bark -     /* is less than .5  */
  509.        ltg[power[i].map].bark < 0.5) {          /* bark, take the   */
  510.        if(power[power[i].next].x > power[i].x ){/* maximum          */
  511.           if(old == STOP) *tone = power[i].next;
  512.           else power[old].next = power[i].next;
  513.           power[i].type = FALSE;
  514.           power[i].x = DBMIN;
  515.           i = power[i].next;
  516.        }
  517.        else {
  518.           power[power[i].next].type = FALSE;
  519.           power[power[i].next].x = DBMIN;
  520.           power[i].next = power[power[i].next].next;
  521.           old = i;
  522.        }
  523.     }
  524.     else {
  525.       old = i;
  526.       i = power[i].next;
  527.     }
  528.  }
  529. }
  530.  
  531. /****************************************************************/
  532. /*
  533. /*        This function calculates the individual threshold and
  534. /* sum with the quiet threshold to find the global threshold.
  535. /*
  536. /****************************************************************/
  537.  
  538. void threshold(power, ltg, tone, noise, bit_rate)
  539. mask FAR power[HAN_SIZE];
  540. g_thres FAR *ltg;
  541. int *tone, *noise, bit_rate;
  542. {
  543.  int k, t;
  544.  double dz, tmps, vf;
  545.  
  546.  for(k=1;k<sub_size;k++){
  547.     ltg[k].x = DBMIN;
  548.     t = *tone;          /* calculate individual masking threshold for */
  549.     while(t != LAST){   /* components in order to find the global     */
  550.        if(ltg[k].bark-ltg[power[t].map].bark >= -3.0 && /*threshold (LTG)*/
  551.           ltg[k].bark-ltg[power[t].map].bark <8.0){
  552.           dz = ltg[k].bark-ltg[power[t].map].bark; /* distance of bark value*/
  553.           tmps = -1.525-0.275*ltg[power[t].map].bark - 4.5 + power[t].x;
  554.              /* masking function for lower & upper slopes */
  555.           if(-3<=dz && dz<-1) vf = 17*(dz+1)-(0.4*power[t].x +6);
  556.           else if(-1<=dz && dz<0) vf = (0.4 *power[t].x + 6) * dz;
  557.           else if(0<=dz && dz<1) vf = (-17*dz);
  558.           else if(1<=dz && dz<8) vf = -(dz-1) * (17-0.15 *power[t].x) - 17;
  559.           tmps += vf;        
  560.           ltg[k].x = add_db(ltg[k].x, tmps);
  561.        }
  562.        t = power[t].next;
  563.     }
  564.  
  565.     t = *noise;        /* calculate individual masking threshold  */
  566.     while(t != LAST){  /* for non-tonal components to find LTG    */
  567.        if(ltg[k].bark-ltg[power[t].map].bark >= -3.0 &&
  568.           ltg[k].bark-ltg[power[t].map].bark <8.0){
  569.           dz = ltg[k].bark-ltg[power[t].map].bark; /* distance of bark value */
  570.           tmps = -1.525-0.175*ltg[power[t].map].bark -0.5 + power[t].x;
  571.              /* masking function for lower & upper slopes */
  572.           if(-3<=dz && dz<-1) vf = 17*(dz+1)-(0.4*power[t].x +6);
  573.           else if(-1<=dz && dz<0) vf = (0.4 *power[t].x + 6) * dz;
  574.           else if(0<=dz && dz<1) vf = (-17*dz);
  575.           else if(1<=dz && dz<8) vf = -(dz-1) * (17-0.15 *power[t].x) - 17;
  576.           tmps += vf;
  577.           ltg[k].x = add_db(ltg[k].x, tmps);
  578.        }
  579.        t = power[t].next;
  580.     }
  581.     if(bit_rate<96)ltg[k].x = add_db(ltg[k].hear, ltg[k].x);
  582.     else ltg[k].x = add_db(ltg[k].hear-12.0, ltg[k].x);
  583.  }
  584. }
  585.  
  586. /****************************************************************/
  587. /*
  588. /*        This function finds the minimum masking threshold and
  589. /* return the value to the encoder.
  590. /*
  591. /****************************************************************/
  592.  
  593. void II_minimum_mask(ltg,ltmin,sblimit)
  594. g_thres FAR *ltg;
  595. double FAR ltmin[SBLIMIT];
  596. int sblimit;
  597. {
  598.  double min;
  599.  int i,j;
  600.  
  601.  j=1;
  602.  for(i=0;i<sblimit;i++)
  603.     if(j>=sub_size-1)                   /* check subband limit, and       */
  604.        ltmin[i] = ltg[sub_size-1].hear; /* calculate the minimum masking  */
  605.     else {                              /* level of LTMIN for each subband*/
  606.        min = ltg[j].x;
  607.        while(ltg[j].line>>4 == i && j < sub_size){
  608.        if(min>ltg[j].x)  min = ltg[j].x;
  609.        j++;
  610.     }
  611.     ltmin[i] = min;
  612.  }
  613. }
  614.  
  615. /*****************************************************************/
  616. /*
  617. /*        This procedure is called in musicin to pick out the
  618. /* smaller of the scalefactor or threshold.
  619. /*
  620. /*****************************************************************/
  621.  
  622. void II_smr(ltmin, spike, scale, sblimit)
  623. double FAR spike[SBLIMIT], scale[SBLIMIT], ltmin[SBLIMIT];
  624. int sblimit;
  625. {
  626.  int i;
  627.  double max;
  628.                 
  629.  for(i=0;i<sblimit;i++){                     /* determine the signal   */
  630.     max = 20 * log10(scale[i] * 32768) - 10; /* level for each subband */
  631.     if(spike[i]>max) max = spike[i];         /* for the maximum scale  */
  632.     max -= ltmin[i];                         /* factors                */
  633.     ltmin[i] = max;
  634.  }
  635. }
  636.         
  637. /****************************************************************/
  638. /*
  639. /*        This procedure calls all the necessary functions to
  640. /* complete the psychoacoustic analysis.
  641. /*
  642. /****************************************************************/
  643.  
  644. void II_Psycho_One(buffer, scale, ltmin, fr_ps)
  645. short FAR buffer[2][1152];
  646. double FAR scale[2][SBLIMIT], ltmin[2][SBLIMIT];
  647. frame_params *fr_ps;
  648. {
  649.  layer *info = fr_ps->header;
  650.  int   stereo = fr_ps->stereo;
  651.  int   sblimit = fr_ps->sblimit;
  652.  int k,i, tone=0, noise=0;
  653.  static char init = 0;
  654.  static int off[2] = {256,256};
  655.  double *sample;
  656.  DSBL *spike;
  657.  static D1408 *fft_buf;
  658.  static mask_ptr FAR power;
  659.  static g_ptr FAR ltg;
  660.  
  661.  sample = (double *) mem_alloc(sizeof(DFFT), "sample");
  662.  spike = (DSBL *) mem_alloc(sizeof(D2SBL), "spike");
  663.      /* call functions for critical boundaries, freq. */
  664.  if(!init){  /* bands, bark values, and mapping */
  665.     fft_buf = (D1408 *) mem_alloc((long) sizeof(D1408) * 2, "fft_buf");
  666.     power = (mask_ptr FAR ) mem_alloc(sizeof(mask) * HAN_SIZE, "power");
  667.     read_cbound(info->lay,info->sampling_frequency);
  668.     read_freq_band(<g,info->lay,info->sampling_frequency);
  669.     make_map(power,ltg);
  670.     for (i=0;i<1408;i++) fft_buf[0][i] = fft_buf[1][i] = 0;
  671.     init = 1;
  672.  }
  673.  for(k=0;k<stereo;k++){  /* check pcm input for 3 blocks of 384 samples */
  674.     for(i=0;i<1152;i++) fft_buf[k][(i+off[k])%1408]= (double)buffer[k][i]/SCALE;
  675.     for(i=0;i<FFT_SIZE;i++) sample[i] = fft_buf[k][(i+1216+off[k])%1408];
  676.     off[k] += 1152;
  677.     off[k] %= 1408;
  678.                             /* call functions for windowing PCM samples,*/
  679.     II_hann_win(sample);    /* location of spectral components in each  */
  680.     for(i=0;i<HAN_SIZE;i++) power[i].x = DBMIN;  /*subband with labeling*/
  681.     II_f_f_t(sample, power);                     /*locate remaining non-*/
  682.     II_pick_max(power, &spike[k][0]);            /*tonal sinusoidals,   */
  683.     II_tonal_label(power, &tone);                /*reduce noise & tonal */
  684.     noise_label(power, &noise, ltg);             /*components, find     */
  685.     subsampling(power, ltg, &tone, &noise);      /*global & minimal     */
  686.     threshold(power, ltg, &tone, &noise,         /*threshold, and sgnl- */
  687.       bitrate[info->lay-1][info->bitrate_index]/stereo); /*to-mask ratio*/
  688.     II_minimum_mask(ltg, <min[k][0], sblimit);
  689.     II_smr(<min[k][0], &spike[k][0], &scale[k][0], sblimit);        
  690.  }
  691.  mem_free((void **) &sample);
  692.  mem_free((void **) &spike);
  693. }
  694.  
  695. /**********************************************************************/
  696. /*
  697. /*        This module implements the psychoacoustic model I for the
  698. /* MPEG encoder layer I. It uses simplified tonal and noise masking
  699. /* threshold analysis to generate SMR for the encoder bit allocation
  700. /* routine.
  701. /*
  702. /**********************************************************************/
  703.  
  704. /****************************************************************/
  705. /*
  706. /*        Fast Fourier transform of the input samples.
  707. /*
  708. /****************************************************************/
  709.  
  710. void I_f_f_t(sample, power)         /* this function calculates */
  711. double FAR sample[FFT_SIZE/2];  /* an FFT analysis for the  */
  712. mask FAR power[HAN_SIZE/2];     /* freq. domain             */
  713. {
  714.  int i,j,k,L,l=0;
  715.  int ip, le, le1;
  716.  double t_r, t_i, u_r, u_i;
  717.  static int M, MM1, init = 0, N;
  718.  double *x_r, *x_i, *energy;
  719.  static int *rev;
  720.  static double *w_r, *w_i;
  721.  
  722.  x_r = (double *) mem_alloc(sizeof(DFFT2), "x_r");
  723.  x_i = (double *) mem_alloc(sizeof(DFFT2), "x_i");
  724.  energy = (double *) mem_alloc(sizeof(DFFT2), "energy");
  725.  for(i=0;i<FFT_SIZE/2;i++) x_r[i] = x_i[i] = energy[i] = 0;
  726.  if(!init){
  727.     rev = (int *) mem_alloc(sizeof(IFFT2), "rev");
  728.     w_r = (double *) mem_alloc(sizeof(D9), "w_r");
  729.     w_i = (double *) mem_alloc(sizeof(D9), "w_i");
  730.     M = 9;
  731.     MM1 = 8;
  732.     N = FFT_SIZE/2;
  733.     for(L=0;L<M;L++){
  734.        le = 1 << (M-L);
  735.        le1 = le >> 1;
  736.        w_r[L] = cos(PI/le1);
  737.        w_i[L] = -sin(PI/le1);
  738.     }
  739.     for(i=0;i<FFT_SIZE/2;rev[i] = l,i++) for(j=0,l=0;j<9;j++){
  740.        k=(i>>j) & 1;
  741.        l |= (k<<8-j);                
  742.     }
  743.     init = 1;
  744.  }
  745.  memcpy( (char *) x_r, (char *) sample, sizeof(double) * FFT_SIZE/2);
  746.  for(L=0;L<MM1;L++){
  747.     le = 1 << (M-L);
  748.     le1 = le >> 1;
  749.     u_r = 1;
  750.     u_i = 0;
  751.     for(j=0;j<le1;j++){
  752.        for(i=j;i<N;i+=le){
  753.           ip = i + le1;
  754.           t_r = x_r[i] + x_r[ip];
  755.           t_i = x_i[i] + x_i[ip];
  756.           x_r[ip] = x_r[i] - x_r[ip];
  757.           x_i[ip] = x_i[i] - x_i[ip];
  758.           x_r[i] = t_r;
  759.           x_i[i] = t_i;
  760.           t_r = x_r[ip];
  761.           x_r[ip] = x_r[ip] * u_r - x_i[ip] * u_i;
  762.           x_i[ip] = x_i[ip] * u_r + t_r * u_i;
  763.        }
  764.        t_r = u_r;
  765.        u_r = u_r * w_r[L] - u_i * w_i[L];
  766.        u_i = u_i * w_r[L] + t_r * w_i[L];
  767.     }
  768.  }
  769.  for(i=0;i<N;i+=2){
  770.     ip = i + 1;
  771.     t_r = x_r[i] + x_r[ip];
  772.     t_i = x_i[i] + x_i[ip];
  773.     x_r[ip] = x_r[i] - x_r[ip];
  774.     x_i[ip] = x_i[i] - x_i[ip];
  775.     x_r[i] = t_r;
  776.     x_i[i] = t_i;
  777.     energy[i] = x_r[i] * x_r[i] + x_i[i] * x_i[i];
  778.  }
  779.  for(i=0;i<FFT_SIZE/2;i++) if(i<rev[i]){
  780.     t_r = energy[i];
  781.     energy[i] = energy[rev[i]];
  782.     energy[rev[i]] = t_r;
  783.  }
  784.  for(i=0;i<HAN_SIZE/2;i++){                     /* calculate power  */
  785.     if(energy[i] < 1E-20) energy[i] = 1E-20;    /* density spectrum */
  786.        power[i].x = 10 * log10(energy[i]) + POWERNORM;
  787.        power[i].next = STOP;
  788.        power[i].type = FALSE;
  789.  }
  790.  mem_free((void **) &x_r);
  791.  mem_free((void **) &x_i);
  792.  mem_free((void **) &energy);
  793. }
  794.  
  795. /****************************************************************/
  796. /*
  797. /*         Window the incoming audio signal.
  798. /*
  799. /****************************************************************/
  800.  
  801. void I_hann_win(sample)             /* this function calculates a  */
  802. double FAR sample[FFT_SIZE/2];  /* Hann window for PCM (input) */
  803. {                                   /* samples for a 512-pt. FFT   */
  804.  register int i;
  805.  register double sqrt_8_over_3;
  806.  static int init = 0;
  807.  static double FAR *window;
  808.  
  809.  if(!init){  /* calculate window function for the Fourier transform */
  810.     window = (double FAR *) mem_alloc(sizeof(DFFT2), "window");
  811.     sqrt_8_over_3 = pow(8.0/3.0, 0.5);
  812.     for(i=0;i<FFT_SIZE/2;i++){
  813.       /* Hann window formula */
  814.       window[i]=sqrt_8_over_3*0.5*(1-cos(2.0*PI*i/(FFT_SIZE/2)))/(FFT_SIZE/2);
  815.     }
  816.     init = 1;
  817.  }
  818.  for(i=0;i<FFT_SIZE/2;i++) sample[i] *= window[i];
  819. }
  820.  
  821. /*******************************************************************/
  822. /*
  823. /*        This function finds the maximum spectral component in each
  824. /* subband and return them to the encoder for time-domain threshold
  825. /* determination.
  826. /*
  827. /*******************************************************************/
  828. #ifndef LONDON
  829. void I_pick_max(power, spike)
  830. double FAR spike[SBLIMIT];
  831. mask FAR power[HAN_SIZE/2];
  832. {
  833.  double max;
  834.  int i,j;
  835.  
  836.  /* calculate the spectral component in each subband */
  837.  for(i=0;i<HAN_SIZE/2;spike[i>>3] = max, i+=8)
  838.     for(j=0, max = DBMIN;j<8;j++) max = (max>power[i+j].x) ? max : power[i+j].x;
  839. }
  840. #else
  841. void I_pick_max(power, spike)
  842. double FAR spike[SBLIMIT];
  843. mask FAR power[HAN_SIZE];
  844. {
  845.  double sum;
  846.  int i,j;
  847.  
  848.  for(i=0;i<HAN_SIZE/2;spike[i>>3] = 10.0*log10(sum), i+=8)
  849.                                                    /* calculate the      */
  850.  for(j=0, sum = pow(10.0,0.1*DBMIN);j<8;j++)       /* sum of spectral   */
  851.    sum += pow(10.0,0.1*power[i+j].x);              /* component in each  */
  852. }                                                  /* subband from bound */
  853. #endif
  854. /****************************************************************/
  855. /*
  856. /*        This function labels the tonal component in the power
  857. /* spectrum.
  858. /*
  859. /****************************************************************/
  860.  
  861. void I_tonal_label(power, tone)     /* this function extracts   */
  862. mask FAR power[HAN_SIZE/2];     /* (tonal) sinusoidals from */
  863. int *tone;                          /* the spectrum             */
  864. {
  865.  int i,j, last = LAST, first, run;
  866.  double max;
  867.  int last_but_one= LAST;
  868.  
  869.  *tone = LAST;
  870.  for(i=2;i<HAN_SIZE/2-6;i++){
  871.     if(power[i].x>power[i-1].x && power[i].x>=power[i+1].x){
  872.        power[i].type = TONE;
  873.        power[i].next = LAST;
  874.        if(last != LAST) power[last].next = i;
  875.        else first = *tone = i;
  876.        last = i;
  877.     }
  878.  }
  879.  last = LAST;
  880.  first = *tone;
  881.  *tone = LAST;
  882.  while(first != LAST){                /* conditions for the tonal     */
  883.     if(first<3 || first>250) run = 0; /* otherwise k+/-j will be out of bounds*/
  884.     else if(first<63) run = 2;        /* components in layer I, which */
  885.     else if(first<127) run = 3;       /* are the boundaries for calc.   */
  886.     else run = 6;                     /* the tonal components          */
  887.     max = power[first].x - 7;
  888.     for(j=2;j<=run;j++)  /* after calc. of tonal components, set to loc.*/
  889.        if(max < power[first-j].x || max < power[first+j].x){   /* max   */
  890.           power[first].type = FALSE;
  891.           break;
  892.        }
  893.     if(power[first].type == TONE){    /* extract tonal components */
  894.        int help=first;
  895.        if(*tone == LAST) *tone = first;
  896.        while((power[help].next!=LAST)&&(power[help].next-first)<=run)
  897.           help=power[help].next;
  898.        help=power[help].next;
  899.        power[first].next=help;
  900.        if((first-last)<=run){
  901.           if(last_but_one != LAST) power[last_but_one].next=first;
  902.        }
  903.        if(first>1 && first<255){     /* calculate the sum of the */
  904.           double tmp;                /* powers of the components */
  905.           tmp = add_db(power[first-1].x, power[first+1].x);
  906.           power[first].x = add_db(power[first].x, tmp);
  907.        }
  908.        for(j=1;j<=run;j++){
  909.           power[first-j].x = power[first+j].x = DBMIN;
  910.           power[first-j].next = power[first+j].next = STOP; /*dpwe: 2nd was .x*/
  911.           power[first-j].type = power[first+j].type = FALSE;
  912.        }
  913.        last_but_one=last;
  914.        last = first;
  915.        first = power[first].next;
  916.     }
  917.     else {
  918.        int ll;
  919.        if(last == LAST) ; /* *tone = power[first].next; dpwe */
  920.        else power[last].next = power[first].next;
  921.        ll = first;
  922.        first = power[first].next;
  923.        power[ll].next = STOP;
  924.     }
  925.  }
  926. }                        
  927.                                 
  928. /****************************************************************/
  929. /*
  930. /*        This function finds the minimum masking threshold and
  931. /* return the value to the encoder.
  932. /*
  933. /****************************************************************/
  934.  
  935. void I_minimum_mask(ltg,ltmin)
  936. g_thres FAR *ltg;
  937. double FAR ltmin[SBLIMIT];
  938. {
  939.  double min;
  940.  int i,j;
  941.  
  942.  j=1;
  943.  for(i=0;i<SBLIMIT;i++)
  944.     if(j>=sub_size-1)                   /* check subband limit, and       */
  945.        ltmin[i] = ltg[sub_size-1].hear; /* calculate the minimum masking  */
  946.     else {                              /* level of LTMIN for each subband*/
  947.        min = ltg[j].x;
  948.        while(ltg[j].line>>3 == i && j < sub_size){
  949.           if (min>ltg[j].x)  min = ltg[j].x;
  950.           j++;
  951.        }
  952.        ltmin[i] = min;
  953.     }
  954. }
  955.  
  956. /*****************************************************************/
  957. /*
  958. /*        This procedure is called in musicin to pick out the
  959. /* smaller of the scalefactor or threshold.
  960. /*
  961. /*****************************************************************/
  962.  
  963. void I_smr(ltmin, spike, scale)
  964. double FAR spike[SBLIMIT], scale[SBLIMIT], ltmin[SBLIMIT];
  965. {
  966.  int i;
  967.  double max;
  968.                 
  969.  for(i=0;i<SBLIMIT;i++){                      /* determine the signal   */
  970.     max = 20 * log10(scale[i] * 32768) - 10;  /* level for each subband */
  971.     if(spike[i]>max) max = spike[i];          /* for the scalefactor    */
  972.     max -= ltmin[i];
  973.     ltmin[i] = max;
  974.  }
  975. }
  976.         
  977. /****************************************************************/
  978. /*
  979. /*        This procedure calls all the necessary functions to
  980. /* complete the psychoacoustic analysis.
  981. /*
  982. /****************************************************************/
  983.  
  984. void I_Psycho_One(buffer, scale, ltmin, fr_ps)
  985. short FAR buffer[2][1152];
  986. double FAR scale[2][SBLIMIT], ltmin[2][SBLIMIT];
  987. frame_params *fr_ps;
  988. {
  989.  int stereo = fr_ps->stereo;
  990.  the_layer info = fr_ps->header;
  991.  int k,i, tone=0, noise=0;
  992.  static char init = 0;
  993.  static int off[2] = {256,256};
  994.  double *sample;
  995.  DSBL *spike;
  996.  static D640 *fft_buf;
  997.  static mask_ptr FAR power;
  998.  static g_ptr FAR ltg;
  999.  
  1000.  sample = (double *) mem_alloc(sizeof(DFFT2), "sample");
  1001.  spike = (DSBL *) mem_alloc(sizeof(D2SBL), "spike");
  1002.             /* call functions for critical boundaries, freq. */
  1003.  if(!init){ /* bands, bark values, and mapping              */
  1004.     fft_buf = (D640 *) mem_alloc(sizeof(D640) * 2, "fft_buf");
  1005.     power = (mask_ptr FAR ) mem_alloc(sizeof(mask) * HAN_SIZE/2, "power");
  1006.     read_cbound(info->lay,info->sampling_frequency);
  1007.     read_freq_band(<g,info->lay,info->sampling_frequency);
  1008.     make_map(power,ltg);
  1009.     for(i=0;i<640;i++) fft_buf[0][i] = fft_buf[1][i] = 0;
  1010.     init = 1;
  1011.  }
  1012.  for(k=0;k<stereo;k++){    /* check PCM input for a block of */
  1013.     for(i=0;i<384;i++)     /* 384 samples for a 512-pt. FFT  */
  1014.        fft_buf[k][(i+off[k])%640]= (double) buffer[k][i]/SCALE;
  1015.     for(i=0;i<FFT_SIZE/2;i++)
  1016.        sample[i] = fft_buf[k][(i+448+off[k])%640];
  1017.     off[k] += 384;
  1018.     off[k] %= 640;
  1019.                         /* call functions for windowing PCM samples,   */
  1020.     I_hann_win(sample); /* location of spectral components in each     */
  1021.     for(i=0;i<HAN_SIZE/2;i++) power[i].x = DBMIN;   /* subband with    */
  1022.     I_f_f_t(sample, power);              /* labeling, locate remaining */
  1023.     I_pick_max(power, &spike[k][0]);     /* non-tonal sinusoidals,     */
  1024.     I_tonal_label(power, &tone);         /* reduce noise & tonal com., */
  1025.     noise_label(power, &noise, ltg);     /* find global & minimal      */
  1026.     subsampling(power, ltg, &tone, &noise);  /* threshold, and sgnl-   */
  1027.     threshold(power, ltg, &tone, &noise,     /* to-mask ratio          */
  1028.       bitrate[info->lay-1][info->bitrate_index]/stereo);
  1029.     I_minimum_mask(ltg, <min[k][0]);
  1030.     I_smr(<min[k][0], &spike[k][0], &scale[k][0]);        
  1031.  }
  1032.  mem_free((void **) &sample);
  1033.  mem_free((void **) &spike);
  1034. }
  1035.